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ABSTRACT 

We present a numerical technique to compute the gravitational lensing induced by sim- 
ulated haloes. It relies on a 2D-Tree domain decomposition in the lens plane combined 
with a description of N-Body particles as extended clouds with a non-singular density. 
This technique is made fully adaptive by the use of a density-dependent smoothing 
which allows one to probe the lensing properties of haloes from the densest regions 
in the center or in substructures to the low-density regions in the outskirts. 'Smooth 
Particle Lensing' has some promising features. First, the deflection potential, the de- 
flection angles, the convergence and the shear are direct and separate end-products 
of the SPL calculation and can be computed at an arbitrary distribution of points 
on the lens plane. Second, this flexibility avoids the use of interpolation or a finite 
differentiation procedure on a grid, does not require padding the region with zeros 
and focuses the computing power on relevant regions. The SPL algorithm is tested by 
populating isothermal spheres and ellipsoids with particles and then comparing the 
lensing calculations to the classical FFT-based technique and analytic solutions. We 
assess issues related to the resolution of the lensing code and the limitations set by the 
simulations themselves. We conclude by discussing how SPL can be used to predict 
the impact of substructures on strong lensing and how it can be generalized to weak 
lensing and cosmic shear simulations. 
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1 INTRODUCTION 



of 



cosmology has now been 

agreeme nt with an exten- 

Tobes ('W ang et al.l |2000|; 

Bahcall et al.l Il999l : 

Spergel et al.l 120031) . One of the central ingredient to 



The concordance model 
widely accepted as it is in 
sive range of observat ional 
Van Waerbeke fc MellieJ 



the model is cold dark matter (CDM). However, in order 
to draw detailed comparison between theoretical models of 
CDM and the observed Universe high resolution computa- 
tional simulations are required. The simplest simulations 
calculate the evolution of density perturbations in a Uni- 
verse, where the matter contribution is composed entirely 
of CD M particl es (Bode & Ostriker 2003; Moore et al. 
p2^ [Navarro et alj [l^; JEvrard <^ Yirg;o Collaboration! 
Il999h . The density distributions that these simulations 
produce cover a wide range of scales, with an enormous 
wealth of information. The challenge that we favirialce now 
is finding ways to efficiently mine this data and compare it 
to observational probes. 

Compact collapsed objects, such as galaxies, groups and 
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clusters, have been studied extensively in the simulations. In 
doing so, a discrepancy k nown as 'the substruct ure prob- 
lem' has been ide ntified ( Moore et al.l [l999l : iKlvpin et al.l 
I1999I : IWeller et al.li2QQ4 '). where the number of small scale 
objects formed in simulated Milky Way's haloes is signifi- 
cantly greater than the number of observed satellites. There 
is some ambiguity as to how dwarf galaxies are matched 
with subhalos of a given mass (see e.g. IStoehr et al.l (|2002l) . 
iKazantzidis et al.l (|200J)). It is perhaps possible to resolve 
the problem for some mass range by adjusting this relation, 
but it is always the case that the simulations predict more 
subhalos than there are observed dwarf galaxies. Whether 
this is a result of suppressed star formation in small halos 
or a deeper problem with the CDM model is not yet clear. 

The flux ratios between images in strong gravita- 
tional lenses are powerful tools for studying substructure 
even wh en it is not assoc i ated w ith any observable stars 
or gas (iMao fc SchneTde3 Il998l: iMetcalf fc MadaiJ | 200i| : 



'Metcalf fc Zhao"2002': 'Dalai fc Kochanek"2nn2': IChibal l2G 
Moustakas fc Metcalf 2003 ; Metcalf et al. 2004- A number 
of authors have investigated the implications of N-Body re- 
sults on strong lensing. One approach is to divide the study 
into two steps. The first step is to study the statistical prop- 
erties of N-Body halo, such as the radial distribution of sub- 
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structure. The second step is to predict the lensing i mphca- 
tions separately by using analytic calculations (Ma o et alJ 
I200A . This method has a number of drawbacks. For in- 
stance, they are dependent on the selection criteria used to 
identify substructure and analytic models for substructure 
lensing do not incorporate all the necessary effects. Another 
approach has been to study the lensing properties of the 
simulations directly using ray-shooting. This approach also 
suffers from a number of drawbacks, which we will investi- 
gate here. 

Both of these methods are limited by the resolutions 
of the N-Body simulations. First, there is a mass resolu- 
tion, set by the fact that halos are modeled using a finite 
number of particles. The current generation of simulations 
typically provide dark matter halos with 1^10 million par- 
ticles within the virial radius. For a Milky Way-sized galaxy, 
this corresponds to mass resolution of 10^ -^ 10^ Mq. Any 
results below this mass scale can only be reached by ex tend- 
ing what we see for larger mass scales (|Mao et al.ll2004ri . Nu- 
merical simulations are also limited in spatial length scales. 
The current generation of simulations have a spatial resolu- 
tion of 50 -^ 500 pc for a galaxy sized halo. Although these 
scales result from an astonishing dynamic range for cosmo- 
logical simulations, they still place major constraints on the 
strong lensing studies since strong lensing probes the very 
inner regions of galaxies, typically radii of 1-10 kpc, . 

Probing these inner regions in detail is also limited by 
the fact that the baryonic component of the galaxy plays 
an important role. An ideal simulation would included dark 
matter as well as the gas, stars and their complex inter- 
actions. Strong l ensing studies of s uch simulations have 
been performed (JBradac et al.ll2004l) . However, it is well- 
known that the current generation of simulations that con- 
tain baryons struggle to reproduce realistic galaxies (see e.g. 
iNav arro et al. (1995), Navarro & Steinmetz (1997)). An al- 
ternative approach is to take the results of an N-Body dark 
matter simulation and add a model for a realistic galaxy 
and then allow the dark matter to relax in the new poten- 
tial ("Amara et all (120 06^). Although somewhat ad hoc, this 
approach allows us to decouple the contributions due to ob- 
served baryonic distributions from those due to the dark 
matter in their halos. 

A final limitation of directly studying the strong 
lensing of N-body halos is the technique used to perform 
the ray shooting. A popular method for calculating the 
lensing properties of object is to project it onto a 2D grid 
and to calculate the gravitational potential using an Fast 
Fourier Transform (FFT) technique. This method has the 
advantage of being significantly faster to compute than 
others, scaling as NlogN and not N^ as in the case of 
particle-particle calculations. The FFT method, however, 
does suffer from two major drawbacks. The first is that 
the resolution of the lensing calculation is limited to grid 
size on which the FFT is performed, and the second is 
that using FFT routines introduces repeating boundary 
conditions into the problem. The latter means that the 
galaxy being studied is no longer an isolated object (see 
Fig. 0. It becomes a part of an endless 2D lattice of 
identical galaxies. The impact of this lattice of external 
galaxies can be controlled by adding a buffer of zeros 
around the galaxy being studied. However, adding a large 
buffer compounds the resolution problem, since for a fixed 




Figure 1. Comparison of the deflection angles calculated using 
the FFT method (blue triangles) and SPL method (red diamonds) 
with the analytic prediction (black curve). On the left (panel (a)) 
we see that for large radii the FFT method becomes dominated 
by the repeating boundary conditions. This is a problem that 
does not occur using the SPL method. However, we know that 
the strong lensing region is very close to the centre. On the left 
(panel(b)) we see that the problem of repeatied boundary condi- 
tions is less pronounced. Here a buffering of one quarter, i.e. one 
quarter of the FFT length is filled with zeros. 



grid size {N x N) the number of grid points that cover the 
galaxy are reduced. Finding the right amount of buffering 
that will balance between these two opposing requirements 
is not a straightforward exercise. Even when the right 
amount of buffering is found, evaluating the impact of these 
constraints on lensing results is not trivial. 

In this paper, we discuss an alternative method for per- 
forming the lensing calculations needed for ray-shooting. We 
focus on the methodology and postpone the applications 
to a forthcoming paper. This technique combines a smooth 
description of particles with well-defined lensing properties 
and a tree based domain decomposition. First we describe 
the principle of ray shooting through FFT and the so-called 
'Smooth Particle Lensing' (SPL hereafter) technique. Then, 
resolution effects and the capacities of SPL are investi- 
gated using realisations of analytic models such as softened 
isothermal spheres and ellipsoids. Finally we conclude by 
discussing the intrinsic limitations induced by simulations 
regarding lensing predictions and present the future appli- 
cations of the SPL technique. 



2 RAY SHOOTING 

To calculate the lensing properties of an object, we first 
adopt the thin lens approximation in which the mass dis- 
tribution is collapsed to a plane perpendicular to the line of 
sight, the lens plane. Light is then assumed to travel along 
straight lines between the source, lens plane and the ob- 
server. Rays experience a discrete deflection at the lensing 
plane. The deflection angle can be calculated directly from 
the mass distribution. This is an excellent approximation for 
galaxy and galaxy cluster size lenses. 
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2.1 Lensing Formalism 

The deflection angle of a lens can be calculated from deflec- 
tion potential, ip{0), of the mass distribution. 



a(0) = \/^{0). 



(1) 



This deflection potential is simply twice the 2D Newtonian 
surface potential and can be calculated from the conver- 
gence, k{0). 



VVW = 2^(6'). 



(2) 



The convergence is dimensionless quantity that can be cal- 
culated from the surface mass distribution on the lensing 
plane, S(^), and geometric factors, 

<o) = £^'-^m = ^, (3) 

where Dl is the angular diameter distance from the observer 
to the lens, Dls is the angular diameter distance from the 
lens to the source, and Ds is the angular diameter distance 
between the observer and the source. These factors are usu- 
ally collected together into a term know as the critical den- 
sity, Sc = c^Ds/4:7tGDlDls- When the surface density of 
a lens is greater than the critical density, ie Ai: > 1, a single 
source can have multiple images. Mapping a position on the 
lensing plane, 0, to the position of the source, /3, is done 
through the lens equation. 



p{0) =0- a{0) 



(4) 



The other properties of a gravitationally lensed system, mag- 
nification and shear, can be calculated using the distortion 
matrix, Aij = —■. Magnification is then /j.{0) = |^|~^5 
which is 



KO) 



■K(e))2-|7(0)|2 



(5) 



(^1 



and the two components of shear are defined as 71 
^22)72 and 72 = A12 = A21. 

In the following sections, we first describe the classical 
method based on the Fourier transforms to compute the 
lensing potential of a simulated halos. Then, we introduce 
the 'Smooth Particle Lensing' technique. 



2.2 FFT 

Solving the Poison equation, (2), can be performed in 
Fourier space, where the equation becomes i'^2p(i) = i^{^), 
where i stands for the Fourier wavenumber. This allows us to 
use Fast Fourier Transforms which use a grid based method 
for finding the Fourier transform of a field. The first step 
is to put the particles onto a 2D grid. Here we do this us- 
ing a Cloud In Cell (CIC) routine. At this stage the first 
undesirable effect of the FFT method is introduced. Plac- 
ing the particles on a grid using CIC introduces a global 
smoothing with an effective CIC window function and lim- 
its the points where the lensing properties are calculated. 
The CIC smoothing is a nuisance in strong lensing because 
this scale does not necessary correspond to scale that is sig- 
nificant to the investigation and so makes the interpretation 
of the results more difficult. It should also be noted that us- 
ing the FFT method for ray shooting is also widely used in 
weak lensing and cosmic shear simulations. In these regimes 



the CIC smoothing is also problematic. When considering 
the weak lensing regime one is typically interested the auto- 
correlations of the lensing field. For such an analysis great 
care must be taken since it is widely known that using a 
simple FFT method to solve the Poisson can result in a loss 
of power on small scales due the CIC smoothing. This loss 
of power affects a large range scales since the CIC window 
function extended in Fourie r space. T his effect is well know 
and can be correct for (Sm ith et alJ 2003 ) , however it is 
not clear in the literature if this correction is applied. In the 
strong lensing regime where the correlation function is usu- 
ally not under investigation this loss of power may be less 
problematic, but the global smoothing is still, never the less, 
undesirable. More problematic in the strong lensing case is 
the fact the lensing calculations are performed on a grid. 

When studying the strong lensing properties of an iso- 
lated halo we are particularly interested in the properties 
of the mass distribution in the very central regions of the 
halo since the impact parameters for light-rays are typically 
small. At the same time, it is important to account for the 
distribution of mass outside this central region since it will 
also have an impact on the lensing signal. This means that 
the entire galaxy must be placed onto a grid and the lensing 
properties calculated at each grid point. However, however 
while locating and analysing images only only the central 
^ 0.1% of the grid points are ever used. This means that 
although the FFT is the fastest known method for com- 
puting the lensing properties at the N grid points, the fact 
that most of these points are not used becomes a signifi- 
cant loss of efficiency. Restricting calculation results to grid 
points is also undesirable, since in order to study multiply 
imaged systems interpolation between the grid points can 
not be avoided since it is highly unlikely that all all images 
will sit on grid points. On small scales, idea lensing solver 
would therefore be able to evaluate the lensing properties of 
a mass distribution at an arbitrary point. For N-Body sim- 
ulations the underlying density distribution is represented 
by point masses. Without any smoothing, the shot noise 
from these point sources would dominate over the signals 
we are studying. Therefore, although a one size fits all CIC 
smoothing is not desirable, a fiexible smoothing is central to 
a lensing solving method. There are two relevant smoothing 
scales. The first is the smoothing scale required to reduce 
shot-noise, which we expect to be linked to the local num- 
ber density of particles. The second smoothing scale is one 
that corresponds to the force resolution of the original N- 
Body simulation. This is because density variation on scales 
smaller than this cannot be trusted to be of real significance. 

On large scales the FFT method also suffers draw backs. 
The introduction of repeating boundary condition has al- 
ready been discussed. The adverse effect of this for strong 
lensing studies is that a buffer needs to be placed around 
the galaxy halo. Since the resolution limit is set by the max- 
imum array that can be manipulated by a CPU, loading 
this array with zeros has a detrimental effect on resolution. 
Once again, the impact of introducing repeating boundary 
conditions on weak lensing and cosmic shear needs to be 
considered carefully since this will impact correlations on 
large scales. It is possible to construct a grid based method 
that does not introduce repeating boundary condition, how- 
ever it complicates the calculation and these methods are 



4 Auhert, Amara & Metcalf 




Figure 2. The properties of each of the particles in the simula- 
tion. Each particle had a 2D Gaussian surface density distribu- 
tion. Top left shows the 2D gravitational deflection angle (a), top 
right shows the convergence (k) and bottom left and right show 
the two components of shear (71 and 72). 



not implemented in the current generation of ray shooting 
routines. 



2.3 Smooth Particle Lensing (SPL) 

Not all lensing calculations in numerical simulations use 
Fourier transform s , for example the method used by 
fWambseanss et al.' (T99oj involves a tree based method to 
calculate deflection angles. Here we suggest to go beyond 
such a calculation in order to compute all the lens quanti- 
ties while using an adaptive description of particles. 

The 'Smooth Particle Lensing' technique is based on the 
on the discrete description of haloes in N-Body simulation 
and relies on the fact that each particle will contribute to the 
projected potential measured at a given point. If i labels the 
i-th particle of a simulated halo, the 2D deflection potential 
measured at r is given by: 



{r) = ^Mr), 



(6) 



where (l)i{r) is the potential created by a single particle. Any 
linear function of the potential can also be expressed as a 
sum over particles. The deflection angle is written as: 



a(r) = V0(r) = ^ V0i(r) = ^ ai(r), 



(7) 



where ai(r) is the 2D deflection angles produced by a single 
particle. Likewise the convergence k and the shear 7 are 



^W = ^Ati(r), 



and 



7(r) = ^7i(r) 



(8) 



(9) 



Therefore, it is possible to derive the lensing properties of an 
arbitrary distribution of points by summing the contribution 
of individual particles. However we are left with two difficul- 
ties at this stage. First, we have to find a convenient single- 
particle lens model, with a 'clean' behaviour from both a 
physical and numerical point of view. Second, we have to 
find an efficient way to perform the sums described in eqs. 
ini Q El s^'iid El These two difficulties are assessed in the next 
two section. 



2.3.1 Smoothed Particles 

A natural way to describe the particles would be to describe 
their density as Dirac's 5{r) functions, leading to 2D poten- 
tial (/)(r) ~ log(r) and deflection angles \a{r)\ ^ 1/r. Not 
only such a choice would induce a noisy projected density 
but it creates singularities in the deflection angles, when 
close-encounters between a ray and a particle occur. An- 
other choice would be to smooth the force appUed to rays 
by the particles, e.g. by assuming \a{r)\ ^ 1/y/r^ + e^. It 
would correspond to the procedure applied in N-Body cal- 
culations. However, it still leads to singular convergence and 
shear for a single particle, making the method less power- 
ful than it can be. Hence, we decided to start from a given 
functional form of a single particle density profile (i.e. its 
convergence /€), to construct its potential back by solving 
the Poisson equation, leading to expressions for the single- 
particle deflection angles and shear. 

Let us consider a particle with a mass rup and a lens 
configuration such as the deflector's redshift is zl and the 
corresponding critical density is given by He. Positions on 
the lens plane are given in terms of physical radii r while 
angles can be recovered by 

'=D^y ^''^ 

where Dl{zl) stands for the angular- diameter distance of 
the lens. We chose to model the single-particle projected 
density by a 2D isotropic 'Gaussian' function (see also figure 
I2J. The related convergence is given by 



K{r) = 



m„ 



27rcr2S< 



■ exp 



r 
'2^ 



(11) 



where r stands for the distance between the ray and the 
particle. The corresponding potential is given by 



(r) 



4.s;(>°s(4^^ 



^^'(-2^))' 



(12) 



where Ei(x) 
given by 

a{r) — rrip — 



TvrT^c 



J°° exp(— x)/xdx. The deflection angle is 



(13) 



One can see in figure that the deflection angles (or equiv- 
alently the force applied to the ray by the particle) has a 
'hollow' behaviour because of the finite but non-nil exten- 
sion of the particle: the force rises as the ray gets closer to 
the 'border' (typically r ^ 3a) of the particle then it adopts 
a behaviour in 1/r as expected. The associated shear created 
by a single particle is given by: 



7i{x,y) = K{r) 



{x'-y') 



r(r,cr) 



(14) 
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and 

j2{x,y) = 2K{r)^r{r,a), 

where 

r(r,cr) =r^ + 2 il-ef^ ] a^ 



center of mass and P satisfies rp < Vo 



If the criterion 



(15) 



(16) 



The a parameter acts as a smoothing parameter and as 
cr ^ 0, the single particle density profile tends to a Dirac- 
5{r) function. This smoothing parameter can be a simple 
constant over all the particles, or in a more sophisticated 
way a function of the considered particle or ray. We show in 
section 13 how this choice can affect the final results. 



2.3.2 Tree 2D 

The single-particle model being set, an efficient way to 
perform the sums over all the particles' contributions re- 
mains to be found. Evidently it cannot be performed by 
direct summation, since the CPU consumption would scale 
as Arrays X A^part • Furthermore an improvement in the sim- 
ulation resolution, i.e. in A^part, would strongly affect the 
computation time, while we want the simulation's resolu- 
tion to be an advantage and not a limitation. For these 
reasons, we suggest to use a summation technique based 
on a tree based domain decomposition. These techniques 
have been widely used in numerical simulations to com- 



pute the 3D f orce created by an a r bitrary distribution of 
Barnes k " """"^ 



particles (e.g, 
0988.) 



3d by an a i 
H^fT986), 
Dikaiakos fc StadeFflQQaV 



Bouchet & Hernauist 
Bertschinger (1998) V 



They are know to scale in a logarithmic way, i.e. they scale 
as Arrays X log A^part in the current case. Here we apply the 2D 
version of this algorithm to compute all the lensing quanti- 
ties created by a set of particles distributed within the lens 
plane. Most of the follo wing has been strongly insp ired by 
the description made bv lDikaiakos fc StadeJ (|l996i) for the 
N-Body code PKDGrav. 

The particles' distribution are organi sed foll o wing a 2D 
version of KD-Tree structures (see e.g. iMoord ([l99l|)) or 
2D Tree. First, the space is divided in two subregions that 
contain the same number of particles. The whole box (or 'the 
root') is connected to two 'branches', which correspond to 
these two subregions. The same 'space-splitting' procedure 
is then applied to these smaller regions (or 'cells'), adding a 
new level to the tree with four new branches. In the end, the 
recursive application of this procedure leads to a binary tree, 
with final branches (or 'leaves') containing A^ieaves particles. 
We chose Meaves = 8. While building the tree, each cell is 
being assigned an opening radius given by 

2rc 



6*^3 



+ rc 



(17) 



The quantities rcom and reenter stand for the maximal pos- 
sible distance between a particle inside the cell and respec- 
tively its center of mass and its geometrical center. The free 
parameter controls the opening radius, as explained below. 
The lensing computation is performed by walking down 
the tree. We discuss this walking procedure for the computa- 
tion of the force felt by a ray (or equivalently the deflection 
angles a) but it translates to the k and 7 computation. Let 
us consider the force felt by a ray at the point P on the lens 
plane. A given cell is opened if the distance rp between its 



is satisfied, the same test is performed on the two branches 
connected to the current cell. This recursion can be stopped 
in two manners. First, the current cell is a leaf : the force 
felt by P is obtained by adding the forces created by each 



particles of this cell. Second possibility, rp > Vo 



the 



force felt by P is obtained by computing the force applied 
by the whole cell, assuming a monopole with a mass equal 
to the number of particles inside the cell and centered on 
its center of mass. In other words, the effect of distant cells 
is modelled as the effect of a macro-particle with the same 
properties as a single one but scaled to the correct mass. 

The parameter appears as a performance control pa- 
rameter. In Eq. 1171 a small implies that all the cells are 
likely to be opened. In this case, the force computation is 
close to the direct summation of the interactions created by 
all the particles, resulting in slow computations. Conversely, 
a large would speed up the procedure but would also lower 
the accuracy of the computation. 

At this stage, we end up with two free parameters: 
which controls the overall summation performance and a 
which introduce a finite spatial resolution. In the following, 
we choose ^ = 0.7: it led to accurate results while ensuring a 
good performance of the summation. The influence of sigma 
is more complex to establish and at the heart of the lensing 
computation at high resolution. 



2.3.3 Adaptive Smoothing 

The final step to the SPL's setting is the implementation of 
an adaptive smoothing. It is clear that if particles sample 
the projected density, their 'extension' should be density- 
dependent. A constant smoothing length may result in over- 
smoothing in high density regions, while it would induce 
shooting noise in low density environment (see also sec. l3.1|) . 
Hence, a should be a function of the local density of par- 
ticles. Two strategies are possible, the first being that each 
dark matter particle is assigned a a by computing its lo- 
cal density. In practice, this involves the computation of the 
local density for 10^ — 10^ particles, which can be highly 
CPU consuming. The second strategy consists of assigning 
the same a to all the particles while its value depend on the 
local density of particles at the position of the light ray. For 
different light- rays a changes. This physically makes sense 
since the density should be correctly estimated at the lo- 
cation of the light-rays, not at the dark matter particles 
location. Because the number of light rays is much smaller 
than the typical number of particles (at least an order of 
magnitude) , this strategy significantly improves the compu- 
tational efficiency. It naturally implies that particles distant 
to the current ray have inappropriate a. However, they only 
contribute to the deflection angle and, being distant, their 
influence varies as rup/r independently of a, no loss of accu- 
racy can be detected. We decide to adopt this second strat- 
egy based on the density at light-rays. 

In our implementation, we chose to estimate this den- 
sity by computing the distance rgph to the A^sph-th closest 
particle to the ray. In practice, this computation is made 
while walking the tree during the lensing computation. The 
density is simply estimated by : 
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iVsph 



(18) 



sph 



This procedure is equivalent to an SPH-density estimation 
with a top-hat Kerne l. Even though more complex kernels 
may be used (see e.g. iLi et alJ (|200(tl ). this simple choice is 
sufficient. Finally, we set a smoothing length for a given ray 
by applying the following relation: 



N^ 



TTyOray 



(19) 



The Na parameter controls the smoothing strength and ap- 
pears as a number of particle over which the smoothing is 
applied. Typically, we found that N^ ~ 64 — 256 give good 
results. 



2.3.4 S PL features 

Having described the technical aspects of SPL, we discuss 
briefly the interesting features of the method. First, a whole 
set of lensing quantities {ip, a, k and 7) are directly com- 
puted from the particles' distribution and this can be ex- 
tended to any Unear function of the potential. Second, a key 
aspect of SPL is its flexibility in terms of geometry. No grid 
or any a priori geometry is required for the sampling of the 
projected density. This methods takes fully advantage of the 
sampling performed by the simulation itself. At any point 
on the lens plane, we simply sum up all the particles con- 
tribution to any lensing quantity. This eliminates the use of 
interpolation procedures, which do not provide additional 
information compared to a coarser sampling. Furthermore, 
the computing power can be focused on an arbitrary region 
with any geometry. The result is an insensitivity of SPL 
to periodic artefacts, preventing the loss of computational 
power and memory in zero-padded regions (see also Fig.QJ. 
For instance, caustics can be mapped back by computing 
the deflections angle only along the critical curves. Also im- 
age distortions can be investigate at very high resolution by 
focusing the SPL calculation around an image spotted at 
lower resolution. 

These features arise from the decoupling in SPL be- 
tween the lens sampling and the computation of the lens' 
effect on the light rays. In N-Body calculations, the relevant 
quantity is the force exerted by particles on particles. There- 
fore, it makes sense to compute the forces at points where 
the density is sampled. In gravitational lensing, the relevant 
calculation is the influence of particles on light-rays. Hence, 
there is no reason to map the lensing quantities on a reg- 
ular grid or at the particle positions. Grid-based methods 
do not achieve this decoupling easily, while it is natural for 
tree-based methods. This makes tree-based (or more gen- 
erally summation-based) methods much more appropriate 
to perform the computations, compared to N-Body calcula- 
tions, where the pros and cons of grid-based and tree-based 
techniques can still be debated. 



3 TESTS USING ANALYTIC DISTRIBUTIONS 

3.1 Non-Singular Isothermal Sphere 

3.1.1 Model 

To test the impacts finite resolution it is important to com- 
pare our measured lensing properties with analytic expec- 
tations. For our initial test we focus on radially symmetric 
case. We begin by defining a convergence field from a Non- 
Singular Isothermal Sphere (NSIS), 



Kid) 



ve^l^t 



(20) 



Through out this paper we scale this distribution to be con- 
sistent with a halo of mass, Mh = 1 x lO^^M©, within a 
radius, rn = 0.5Mpc. The halo is then assumed to be at a 
red-shift of z=0.8 and the background sources are assumed 
to be at a red-shift of z=3.0. 

Solving the Poisson equation (0) for the convergence 
field described earlier leads to the lensing potential. 



^{0) = 20ono 



'1 + 




and a deflection field of, 



a{e) =2ko-^ 



1 + ^ 



(21) 



(22) 



The subsequent shear and magnification are simply calcu- 
lated from this deflection angle by constructing the distor- 
tion matrix. Ay and using the definition above equation |51 
In order to draw conclusions about current numerical 
simulation we chose to study a NSIS were the core radius is 
2000 times smaller than a typical virial radius for a galaxy 
of mass M. In our case the M is lO^^M©. This core radius 
was chosen because it is approximately the smoothing scale 
expected in current N-body simulations. These simulations 
contain roughly 10^ particles. To accurately simulate an 
NSIS of this type under these conditions, we ffist converted 
the density distribution into a probability distribution func- 
tion (PDF). We then randomly sampled this distribution 
10^ times to obtain positions for the particles. Finally by 
setting the mass of each of these points to 10^ M© we were 
able to produce realisation of our NSIS with particle sam- 
pling similar to that of an N-body halo. 



3.1.2 Radial Profiles 

As an illustration, we begin by testing the impact of smooth- 
ing scale with a constant smoothing scale a. We do this by 
calculating the inverse magnification for a number of rays 
which are placed on concentric rings around the centre of 
the NSIS. Such a test is crucial since magnification is the 
quantity that is most sensitive to variations in density, be- 
cause it is a combination of second order derivatives of the 
lensing potential. Furthermore, magnification serves to spot 
the location of critical lines and caustics. The radii of the 
rings are chosen so that they are distributed evenly in log(^) 
which allows us to cover a wide range of scales and hence 
take full advantage of the sampling flexibility of the SPL 
method. 

Figure El shows the inverse magnification results. For 
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Figure 3. The inverse magnification of a NSIS with a core ra- 
dius of 0.037 arcsec (r/j,/2000). The sohd fine shows the analytic 
solutions. The points are the results of SPL calculations on a dis- 
tribution with 10^ particles. In each case a constant smoothing 
scale is used, the triangles (red) are for a smoothing of e = 0.0074 
arcsec, the diamonds (blue) are for e = 0.074 arcsec and the circles 
(green) are for e = 0.74 arcsec. 



each case the smoothing was fixed to the same value over all 
the points calculated. The green circles, blue diamond and 
red triangles are, respectively, for smoothing scales of 100, 
1000, and 1000 times smaller than the galaxy radius. From 
this we clearly see that a variable smoothing scale would 
be desirable. In the outer roughly 2 arcseconds, we see that 
the green point accurately follow the analytic predictions 
with minimal spread, where as the red triangles suffer a 
substantial spread due to shot noise. In the inner regions 
however the green points are incapable of finding the crit- 
ical curves, points where the inverse magnification is zero, 
hence we would not recognise this as a system capable of 
having multiple images . With the blue points we are able 
to accurately find the radial critical curve (~ 1 arcsec), with 
very little spread, but fail to capture the tangential critical 
curve (^ 0.15 arcsec). In these very inner regions the red 
points do better but the spread is still substantial. 

From now on, we consider the case where an adaptive 
smoothing is set following the procedure described in section 
12.3.31 Unless specified otherwise, we set Na = 256. Figure 3] 
shows the deflection angle calculation for the same halo as 
in figure El As expected we see that the deflection angle is 
easily calculated for a wide range of projected radii. Figure 
1^1 we clearly see that the convergence and the shear fields 
are more sensitive to shot noise than the deflection angles, 
but they are still accurately reproduced. Finally these are 
combined to calculate the inverse magnification shown in 
figure El Here both the results from the SPL method and 
the FFT method are shown. We can see from this figure that 
the SPL does offers greater flexibility than the FFT method. 
We are able to sample the field at points of interest to us, in 
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Figure 4. The defiection angle of a NSIS with a core radius 
of 0.037 arcsec (r/j,/2000). The solid black line is the analytic 
solution and the blue diamond symbols are those measured using 
the SPL method with adaptive smoothing. 



this particular case using concentric rings and we are able to 
control the spread thanks to the adaptive smoothing scale. 
Finally, we show in figure |3 the same calculations per- 
formed on the same model but sampled using only 10^ 
particles. We applied an adaptive smoothing with N^ = 
32, 64, 128, 256. The deflection angles are accurately com- 
puted, even though large smoothing values seem to slightly 
affect the calculation close to the halo's center. On the other 
hand, the inverse magnification is clearly more sensitive to 
the lower resolution of this halo, especially close to the cen- 
ter. Large smoothing value lead to a limited spread but the 
inner rise of the inverse magnification is incorrectly repro- 
duced and occurs at larger radii than in the analytic solu- 
tion. Conversely, small smoothing value seem to provide a 
better fit to the analytic magnification profile in these inner- 
most regions, but the spread become very important. Con- 
sequently, constrains on the radial caustic would be highly 
uncertain even though they would be closer on average to 
the prediction. We conclude that such a halo is not suited to 
study the lensing properties of such models via SPL and we 
argue that it is very unlikely that another method would be 
able to compute correctly the lensing signal of the underly- 
ing model. Consequently, the limitation is clearly induced by 
the halo (or the N-Body simulation) and not by an inability 
to compute the signal properly. 



3.1.3 Magnification Maps 

The radial properties of the lensing fields discussed thus far 
give insight into the strengths is each method. It is also very 
important to look at the lensing field in 2D. Figure |H1 shows 
the inverse magnification for four difference cases. Panel (a) 
shows the results of the highest resolution FFT grid calcu- 
lation that we were able to perform (4096 x 4096) . The pixel 
scale on the this zoom in is clearly visible since the image is 
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Figure 5. The solid black line shows the analytic solution for 
the convergence of the NSIS with a core radius of 0.037 arc- 
sec (r/j,/2000) and dashed line shows the magnitude of the shear 
(I7I). Overplotted are the convergence (blue squared) and the ab- 
solute shear (blue crosses) calculated using the SPL method on 
a distribution containing 10^ particles using adaptive smoothing 
(Na = 256). 




Figure 6. The inverse magnification of the NSIS with a radius of 
0.037 arcsec (r/i/2000). The solid black curve shows the analytic 
solution. In red are the results of an FFT calculation performed on 
a 4096"^ grid on a distribution made up of 10^ particles. The blue 
diamonds show the SPL measurements on the same distribution 
using adaptive smoothing {Na = 256). 
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Figure 7. Deflection angle (Top) and inverse magniflcation (Bot- 
tom) of the NSIS with a core radius of 0.037 arcsec (r/,/2000). 
The model is sampled with le6 particles. The SPL calculation is 
performed with N^ = 256, 128, 64, 32 (resp. gold, blue, red and 
green points). The solid lines stand for the analytic proflles. 



made up of ^ 160 x 160 pixels. Although it may be possible 
to optimise our routines further it is not likely that we will 
be able to increase the resolution by more than a factor of 
2. Panel (b) of figure |H1 shows a set of SPL calculations per- 
formed at the same grid points as the FFT results of panel(a) 
and with a fixed smoothing scale that is comparable to the 
CIC smoothing. These two 'low resolution' results are in 
good agreement with each other. However, since the SPL 
calculations where only performed on the points shown, the 
SPL method was faster. Panel (c) shows SPL results where 
the inverse magnification is sampled with 1024 x 1024 points 
on a polar grid with a fixed smoothing scale of 10~^ in units 
of the virial radius. Here the features from panels (a) and 
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Figure 8. The inverse magnification for the NSIS. Panel (a) 
shows the results from the FFT calculation on 4096^. Panel (b) 
shows the result from the SPL calculation with a smoothing scale 
comparable to the CIC smoothing scale of the FFT calculation 
with sampling at the same grid points as the FFT calculation of. 
Panel (c) shows the inverse magnification using the SPL routine 
where where the number of points where the lensing calculation 
is increased to 1024 x 1024 pixels and the smoothing scale is set 
to a constant lO"'^. In this panel, coordinates are given in units 
of the virial raidus. The same sampling is used in panel (d) as 
in panel (c), however in (d) an adaptive smoothing scale that is 
linked to underlying surface density is used. 



(b) can been seen in detail and since they are persistent 
throughout all the calculations we must conclude that the 
'flame' patterns are indeed real features unique to this real- 
isation of the halo. They are there a result of shot noise due 
to finite mass resolution. Panel (d) shows the image of the 
inverse magnification when using adaptive smoothing with 
Na = 256. We see that by doing this we have managed to 
greatly reduce the 'flame' features, although the pattern is 
not perfectly symmetric as would be the case with infinite 
mass resolution. 

3.1.4 Critical and Caustic Curves 

The final remaining features of the NSIS are the critical 
and caustic curves. These are shown in figure El Here the 
solid curves are the caustics and the dashed are the critical 
curves. Shown in red are the curves we measure using our 
SPL method and in black are the analytic solutions for our 
NSIS. Clearly the two results agree very well and provide a 
2D validation of the whole technique. 

The location of SPL's critical curves also permits us to 
investigate the halo's mass resolution issues. In Fig. 1101 we 
show the critical lines computed for two haloes with 10^ and 
10 particles and for two different smoothing [Na — 32, 256). 
Considering the 10^ halo first, increasing Na removes the 
small scales features in the magnification map and reduces 
the spread around the analytic solution the critical curves. In 
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Figure 9. The critical and caustic curves for the NSIS. The 
dashed red lines are the critical curves the solid red lines are 
the corresponding caustics. The black curves are the analytic so- 
lutions and the red line are the results for SPL on a distribution 
with 10^ particles and core radius of 0.037 arcsec (r/j,/2000). 



general, the overall structure of these lines is well reproduced 
in all the cases. For the 10^ particle halo, the same behaviour 
occurs for the outer critical curve: more smoothing gives 
reduced spread and a better agreement with the analytic 
solution. However the inner curve drifts outwards from the 
correct location as N^^ increases, even though the spread is 
reduced. Clearly, no good compromise can be found between 
resolution and noise control for a halo with this lower mass 
resolution. 

We stress again that no grid-based interpolation has 
been used here to compute the deflection angle on the critical 
lines and map them back to caustics. A first computation 
has been performed to spot the critical lines (using the same 
procedure which results in panel (d) in Fig. |H1). Then the 
exact deflection angles were computed along these critical 
lines only. 



3.2 Non-Singular Isothermal Ellipsoid 

We continue our investigation of the SPL technique by 
breaking the spherical symmetry. We computed the mag- 
nification maps and the deflection angles for the Non- 
Si ngular Isothermal Ellip soids (NSIE hereafter) described 
bv lKormann et all (|l994 h. These models are defined by the 
following projected density: 



m 



cVe 



2^/^F+^ 



(23) 



where the scaling Sc has been added to allow us to chose 
halo masses configurations. We have selected models with a 
halo of mass Mh = 1 x 10^^ M© with a main axis' length 
Th — 0.5 Mpc. The lens configuration is the same as the 
one chosen for the NSIS, i.e. sources at z = 3 and a lens 
at z = .8. Aside from the total mass and core radius, an 
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Figure 10. Critical curve locations for our template NSIS sam- 
pled with 10^ (left column) and 10^ (right column) particles. SPL 
calculation are performed with No- = 32, 256 (resp. bottom and 
top row) on a 256x256 logarithmically sampled polar grid. The 
dashed lines stand for the analytically derived locations of the 
critical curves. The solid lines stand for the SPL calculations. 



additional parameter e, the axis ratio, introduces ellipticity 
into the model through. 



= y/WTe^\ 



(24) 



where Ox and 6y are the direction in the x — y plane. We 
generated three 10^ particles' models with e = 0.2, 0.4. and 
0.8 and Fig. ITTl shows the critical and caustic curves loca- 
tion on the lens plane. SPL computations were performed 
with Na — 128 on a 256x256 logarithmic polar grid with 
10~^ < r/rh < 5 • 10~^. Because the spherical symmetry is 
broken, the radial caustic is not degenerate anymore: it al- 
lows to probe the accuracy of the outer critical curve's com- 
putation in a regime where the projected density is lower 
than in the central regions. Clearly the match between SPL 
calculations (shown in red) and the analytic solution (shown 
in black) is good and the lines are almost indistinguishable. 
Realistic dark matter haloes exhibits a certain level of tri- 
axiality which should be easily handled with the current 
technique. In particular, 'naked cusps' (i.e. cusps outside 
the radial caustics) are accurately reproduced, even for the 
most flattened models. Therefore image multiplicity statis- 
tics, the relative numbers of two-, four- images regions and 
naked cusp systems, should be well reproduced by the sim- 
ulations. 



4 SUMMARY k PROSPECTS 

We present a method for computing the gravitational lensing 
induced by simulated gravitational lenses. It solves the Pois- 
son equation using a 2D-Tree decomposition technique com- 
bined with a description of simulation particles as 'extended 
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Figure 11. Critical (dashed) and caustic (plain) curve for the 
three different NSIE models: e = 0.2, 0.4, 0.8 (from top to bot- 
tom). SPL calculations (red) are compared to analytic solu- 
tions (black). NSIE models are sampled with 10^ particles and 
Mh = lO^^M© and a core radius, re = r/,/2000 with r/, = 500 
kpc. SPL calculations were performed with Na = 128 on a 
256x256 polar grid. 
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clouds'. This so-called 'Smooth Particle Lensing' technique 
allows the direct computation of any linear function of the 
deflection potential, such as the deflections angles, conver- 
gence, and shear. We test this method using analytic mod- 
els of Non-singular isothermal spheres (NSIS) and ellipsoids 
(NSIE). The lensing properties are very accurately repro- 
duced if the mass resolution is high enough. For instance, 
the caustic structure including naked cusps are correctly re- 
produced even for highly flattened systems and therefore ac- 
curate predictions of the image multiplicities are expected. 

These tests showed that adaptive smoothing is neces- 
sary if one wants to probe the whole lens plane to sufficient 
accuracy: a small smoothing scale does well in high den- 
sity regions but induces a large scatter in low density envi- 
ronments. Conversely, over-extended particles smooth over 
small-scale features such as the inner cusp and substruc- 
tures. In trying to achieve a compromise between smooth- 
ing and resolution we found that the number of particles in 
simulated haloes can be critical to locating the critical and 
caustic curves. For instance, a NSIS with a 'core' as small 
as rh,/2000 and sampled with 10^ particles cannot be fully 
investigated in our fiducial lens configuration (zi = 0.8 and 
Zs = 3.). This emphasises the need of high resolution haloes 
in order to accurately predict the strong lensing produced 
by these objects. 

Because it combines noise-control with high resolution, 
this tool is particularly suited to studying the effect of sub- 
structures in dark matter halos, topics such as violations 
of the cusp caustic relation between image magnifications 
or predictions of image mu ltiplicities (see e.g Bradac et all 
(|2004 ). lAmara et all (^2006)). Because it does not suffer from 
periodic boundary conditions, SPL is naturally extended to 
weak lensing studies of isolated objects such as galaxy clus- 
ters. In a large simulation box the method will adapt to 
put computational power where it is most needed, improv- 
ing the accuracy of shear correlations over a large range of 
scales (from clusters to degree-scale cosmic shear surveys) 
while reducing CPU time. Since a significant fraction of the 
N-body codes are based on tree structures, one could even 
imagine on-the-fly lensing calculations during the dynamical 
integration. A level of accuracy is achieved in our implemen- 
tation of SPL where the calculation is effectively limited by 
the accuracy of the simulation used to find the mass distri- 
bution rather than the lensing code itself. 
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